Linear regressions on spiracle data with non-parametric confidence intervals for slopes¶
To begin with, we need to import necessary python packages.
In [1]:
import numpy as np
import pandas as pd
import statsmodels.api
import scipy.stats
import bokeh.io
import bokeh.plotting
from Bio import Phylo
import io
import cmdstanpy
import arviz as az
import bebi103
bokeh.io.output_notebook()
We can now read in the data into a dataframe for analyis.
In [2]:
df = pd.read_csv("./20190322_supp_table_2.csv")
We take a look at the format for the data.
In [3]:
df['species_underscore'] = [spec.replace(" ", "_") for spec in df['species']]
df.head()
Out[3]:
For some of this analysis, we will look at the per-species averages for our measurements. To get this, we use a simple aggregate function on the dataframe and take a look at the results.
In [4]:
df_averages = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(np.average)
df_averages['subfamily'] = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(max)['subfamily']
df_averages.head()
Out[4]:
In [5]:
species_per_subfam=df_averages.groupby(['subfamily', 'spiracle'], as_index=False).count().groupby('subfamily').aggregate(max).reset_index()[['subfamily', 'species']]
species_per_subfam.columns = ('subfamily', 'subfam_count')
species_per_subfam
Out[5]:
In [6]:
df_averages = df_averages.merge(species_per_subfam, on='subfamily')
For our plots, we will log transform the data. We will add a column to the dataframe with the log transformed data.
In [7]:
df_averages.head()
Out[7]:
In [8]:
df_averages['log area (mm^2)'] = np.log10(df_averages['area (mm^2)'])
df_averages['log dist'] = np.log10(df_averages['depth (mm)'])
df_averages['log mass (g)'] = np.log10(df_averages['mass (g)'])
df_averages['log area/dist'] = np.log10(df_averages['area (mm^2)']/df_averages['depth (mm)'])
df_averages['log area^2/dist'] = np.log10(df_averages['area (mm^2)']**2/df_averages['depth (mm)'])
df_averages.head()
Out[8]:
In addition to log transforming the species averaged data, we will do the same for the whole data set.
In [9]:
df['log area (mm^2)'] = np.log10(df['area (mm^2)'])
df['log dist'] = np.log10(df['depth (mm)'])
df['log mass (g)'] = np.log10(df['mass (g)'])
df['log area/dist'] = np.log10(df['area (mm^2)']/df['depth (mm)'])
df['log area^2/dist'] = np.log10(df['area (mm^2)']**2/df['depth (mm)'])
df.head()
Out[9]:
In [10]:
tree = Phylo.read(io.StringIO("(((Cyclocephala_borealis:1.0,(Dynastes_hercules:1.0,Trypoxylus_dichotomus:1.0):1.0):1.0,Popilia_japonica:1.0):1.0,(Protaetia_orientalis:1.0,((((Coelorrhina_hornimani:1.0,Mecynorrhina_torquata:1.0):1,Dicronorrhina_derbyana:1):1,Eudicella_euthalia:1.0):1,Goliathus_goliathus:1.0):1.0):1.0):0.0;"),
"newick")
In [11]:
Phylo.draw(tree)
In [12]:
species = df.sort_values('mass (g)')['species_underscore'].unique()
cov = np.zeros(shape=(len(species), len(species)))
comp_matrix = np.zeros(shape=(len(species), len(species))).astype(str)
for i, spec in enumerate(species):
for j, comp in enumerate(species):
if i == j:
cov[i, j] = 6.0
continue
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
comp_matrix[i, j] = spec + ':' + comp
In [13]:
cov
Out[13]:
In [14]:
print(pd.DataFrame(cov))
pd.DataFrame(comp_matrix)
Out[14]:
In [15]:
df_averages.loc[df_averages['spiracle'] == '1']
Out[15]:
With our data in this format, we can build the PGLS model for our regression. To do this, we will write a statistical model for the data. We choose a multivariate normal distribution with covariance matrix given by the above matrix adjusted by the parameter $\lambda$ (as in PGLS), and put a prior on $\lambda$. This gives
\begin{align} f(x_i; a, b) &= ax_i + b \\ b &\sim \text{Norm}(0, 1) \\ a &\sim \text{Norm}(0, 2) \\ \lambda &\sim \text{beta}(1.1, 1.1) \\ y_i &\sim \text{multi_normal}\left(f_\mu(x_i; a, b), f_\sigma(\mathrm{cov}, \lambda)\right) \end{align}Time for a prior predictive check on our model!
In [16]:
sm_prior_pred = cmdstanpy.CmdStanModel(
stan_file="spiracle_pgls_prior_predictive.stan"
)
In [17]:
x = df_averages.sort_values('mass (g)').loc[df_averages['spiracle'] == '1', 'log mass (g)'].values
data = {
"N": len(x),
"x": x,
"a_mu": 0.66,
"a_sig": 0.3,
"b_mu": -1.0,
"b_sig": 1.0,
"lambda_alpha": 1.1,
"lambda_beta": 1.1,
"cov_phylo":cov,
}
samples = sm_prior_pred.sample(data=data, sampling_iters=1000, fixed_param=True)
In [18]:
samples = az.from_cmdstanpy(posterior=samples, prior=samples, prior_predictive=['y'])
In [19]:
bokeh.io.show(
bebi103.viz.predictive_regression(
samples.prior_predictive['y'],
samples_x=x,
percentiles=[30, 60, 90, 99],
x_axis_label='mass',
y_axis_label='area',
)
)
In [20]:
p = bokeh.plotting.figure()
p.circle('log mass (g)', 'log area (mm^2)', source=df_averages.loc[df_averages['spiracle'] == '1'])
bokeh.io.show(p)
In [21]:
size = 1000
lines = [a*x+b+sample for a, b, sample in zip(np.random.normal(0,2, size=size), np.random.normal(0,1,size=size),
np.random.multivariate_normal(x, cov, size=size))]#+np.random.normal(0, 3, size=(100,) ))]
In [ ]:
In [22]:
p = bokeh.plotting.figure()
[p.line(x, y, alpha=0.2) for y in lines]
bokeh.io.show(p)
In [23]:
species = df.sort_values('mass (g)')['species_underscore'].unique()
cov = np.zeros(shape=(len(species), len(species)))
comp_matrix = np.zeros(shape=(len(species), len(species))).astype(str)
for i, spec in enumerate(species):
for j, comp in enumerate(species):
#if i == j:
# cov[i, j] = 6.0
# continue
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
comp_matrix[i, j] = spec + ':' + comp
In [24]:
size = 1000
lines = [a*x+b+sample for a, b, sample in zip(np.random.normal(0,2, size=size), np.random.normal(0,1,size=size),
np.random.multivariate_normal(x, cov, size=size))]#+np.random.normal(0, 3, size=(100,) ))]
In [25]:
p = bokeh.plotting.figure()
[p.line(x, y, alpha=0.2) for y in lines]
bokeh.io.show(p)
In [26]:
cov
Out[26]:
In [27]:
plots = []
sample = np.random.multivariate_normal(x, cov, size=1000)
num = df_averages.sort_values('mass (g)').loc[df_averages.sort_values('mass (g)')['spiracle']=='1', 'subfam_count'].values
for i in range(len(sample[0])):
for j in range(len(sample[0]))[i:-1]:
p = bokeh.plotting.figure(plot_width=200, plot_height=200, x_range=(-5, 5), y_range=(-5, 5), title='comp of ' +str(i)+' and '+str(j)+' | '+str(num[i])+','+str(num[j]))
p.circle(sample[:, i], sample[:, j], alpha=0.1)
plots.append(p)
In [28]:
bokeh.io.show(bokeh.layouts.gridplot(plots, ncols=4))
In [29]:
p = bokeh.plotting.figure()
[p.circle(x, sample, alpha=0.05) for sample in np.random.multivariate_normal(x, cov, size=1000)]
bokeh.io.show(p)
In [ ]:
In [30]:
print(df_averages.sort_values('mass (g)').loc[df_averages.sort_values('mass (g)')['spiracle']=='1', 'subfam_count'].values)
In [31]:
Phylo.draw(tree)
In [32]:
species = df.sort_values('mass (g)')['species_underscore'].unique()
cov = np.zeros(shape=(len(species), len(species)))
comp_matrix = np.zeros(shape=(len(species), len(species))).astype(str)
for i, spec in enumerate(species):
for j, comp in enumerate(species):
if i == j:
cov[i, j] = 6.0
continue
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
comp_matrix[i, j] = spec + ':' + comp
In [33]:
cov*1/1000
Out[33]:
In [68]:
x = df_averages.sort_values('mass (g)').loc[df_averages['spiracle'] == 'S', 'log mass (g)'].values
y = df_averages.sort_values('mass (g)').loc[df_averages['spiracle'] == 'S', 'log area/dist'].values
data = {
"N": len(x),
"x": x,
"y": y,
"cov_phylo":cov,
}
In [80]:
sm = cmdstanpy.CmdStanModel(stan_file='spiracle_pgls_model.stan')
samples = sm.sample(data=data, sampling_iters=1000, chains=4)
samples = az.from_cmdstanpy(posterior=samples, posterior_predictive=["y_ppc"])
In [81]:
bebi103.stan.check_all_diagnostics(samples, var_names=['b', 'a'])
Out[81]:
In [83]:
bokeh.io.show(
bebi103.viz.corner(
samples, pars=["a", "b", "lambda", "sigma"], xtick_label_orientation=np.pi / 4
)
)
In [84]:
p = bokeh.plotting.figure()
norm = np.linspace(-1, 1.5, 200)
bebi103.viz.ecdf(samples.posterior['a'].values.ravel(), color='red', p=p)
p.line(norm, scipy.stats.norm.cdf(norm, 0.33, 0.3))
bokeh.io.show(p)
In [85]:
az.summary(samples, var_names=['b', 'a', 'lambda'])
Out[85]:
In [86]:
az.summary(samples, var_names=['b', 'a'])['mean']['a'] + az.summary(samples, var_names=['b', 'a'])['sd']['a']*1.96
Out[86]:
In [87]:
az.summary(samples, var_names=['b', 'a'])['mean']['a'] - az.summary(samples, var_names=['b', 'a'])['sd']['a']*1.96
Out[87]:
In [88]:
bokeh.io.show(bebi103.viz.ecdf(samples.posterior['a'].values.ravel()-0.333))
In [89]:
np.sum(samples.posterior['a'].values.ravel()-0.333 < 0)/len(samples.posterior['a'].values.ravel())
Out[89]:
In [91]:
y_ppc = samples.posterior_predictive['y_ppc'].stack(
{"sample": ("chain", "draw")}
).transpose("sample", "y_ppc_dim_0")
bokeh.io.show(
bebi103.viz.predictive_regression(
y_ppc,
samples_x=x,
percentiles=[30, 50, 70, 95],
data=np.vstack((x, y)).transpose(),
x_axis_label='ant density, k [ant/cm^2]',
y_axis_label='ant flow, q [ant/cm/s]',
x_range=[x.min()-0.1, x.max()+0.1],
)
)
In [92]:
p = bokeh.plotting.figure(x_axis_label='log mass', y_axis_label='log area')
Y = df.sort_values('mass (g)').loc[df['spiracle'] == '2', 'log area/dist'].values
X = df.sort_values('mass (g)').loc[df['spiracle'] == '2', 'log mass (g)'].values
X = statsmodels.api.add_constant(X)
model = statsmodels.api.OLS(Y,X)
results = model.fit()
results.params
p.circle(df.sort_values('mass (g)').loc[df['spiracle'] == '2', 'log mass (g)'].values,
df.sort_values('mass (g)').loc[df['spiracle'] == '2', 'log area/dist'].values, color='black', size=15, alpha=0.3, fill_color='green')
p.line(np.linspace(x.min(), x.max(), 200),
results.params[1]*np.linspace(x.min(), x.max(), 200) + results.params[0],
line_width=5, line_dash='dashed', color='orange', alpha=0.3)
Y = y
X = x
X = statsmodels.api.add_constant(X)
model = statsmodels.api.OLS(Y,X)
results = model.fit()
results.params
p.line(np.linspace(x.min(), x.max(), 200),
az.summary(samples, var_names=['b', 'a'])['mean']['a']*np.linspace(x.min(), x.max(), 200) + az.summary(samples, var_names=['b', 'a'])['mean']['b'],
line_width=5, color='black')
p.line(np.linspace(x.min(), x.max(), 200),
results.params[1]*np.linspace(x.min(), x.max(), 200) + results.params[0],
line_width=5, line_dash='dashed', color='grey')
p.circle(x, y, color='black', size=15, alpha=0.5, fill_color='red')
bokeh.io.show(p)
In [46]:
p = bokeh.plotting.figure(x_axis_label='log mass', y_axis_label='log area')
p.circle(df.sort_values('mass (g)').loc[df['spiracle'] == '1', 'log mass (g)'].values,
df.sort_values('mass (g)').loc[df['spiracle'] == '1', 'log area (mm^2)'].values, color='black', size=15, alpha=0.8, fill_color='red')
bokeh.io.show(p)
In [ ]:
In [47]:
species = df.sort_values('mass (g)')['species_underscore'].unique()
cov = np.zeros(shape=(len(species), len(species)))
comp_matrix = np.zeros(shape=(len(species), len(species))).astype(str)
for i, spec in enumerate(species):
for j, comp in enumerate(species):
comp_matrix[i, j] = spec + ':' + comp
if i == j:
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
continue
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
cov
Out[47]:
In [48]:
cov=np.linalg.inv(cov)
print(pd.DataFrame(np.round(cov, decimals=3)))
In [49]:
x = df_averages.sort_values('mass (g)').loc[df_averages['spiracle'] == '3', 'log mass (g)'].values
y = df_averages.sort_values('mass (g)').loc[df_averages['spiracle'] == '3', 'log area (mm^2)'].values
data = {
"N": len(x),
"x": x,
"y": y,
"cov_phylo":cov,
}
In [50]:
sm = cmdstanpy.CmdStanModel(stan_file='spiracle_pgls_model.stan')
samples = sm.sample(data=data, sampling_iters=1000, chains=4)
samples = az.from_cmdstanpy(posterior=samples, posterior_predictive=["y_ppc"])
In [51]:
bebi103.stan.check_all_diagnostics(samples, var_names=['b', 'a', 'lambda'])
Out[51]:
In [52]:
bokeh.io.show(
bebi103.viz.corner(
samples, pars=["a", "b", "lambda"], xtick_label_orientation=np.pi / 4
)
)
In [53]:
y_ppc = samples.posterior_predictive['y_ppc'].stack(
{"sample": ("chain", "draw")}
).transpose("sample", "y_ppc_dim_0")
bokeh.io.show(
bebi103.viz.predictive_regression(
y_ppc,
samples_x=x,
percentiles=[30, 50, 70, 99],
data=np.vstack((x, y)).transpose(),
x_axis_label='ant density, k [ant/cm^2]',
y_axis_label='ant flow, q [ant/cm/s]',
x_range=[x.min()-0.1, x.max()+0.1],
)
)
In [54]:
p = bokeh.plotting.figure()
p.line(np.linspace(-1, 1.5, 200), np.linspace(-1, 1.5, 200)*0.7672-1.4551)
p.line(np.linspace(-1, 1.5, 200), np.linspace(-1, 1.5, 200)*0.5887-1.3282, color='red')
p.line(np.linspace(-1, 1.5, 200), np.linspace(-1, 1.5, 200)*az.summary(samples, var_names=['b', 'a', 'lambda'])['mean'][1]
+ az.summary(samples, var_names=['b', 'a', 'lambda'])['mean'][0], color='black')
p.circle(x, y, size=15)
bokeh.io.show(p)
In [55]:
az.summary(samples, var_names=['b', 'a', 'lambda'])
Out[55]:
In [56]:
library(nlme)
cat("((((Homo:0.21,Pongo:0.21):0.28,",
"Macaca:0.49):0.13,Ateles:0.62):0.38,Galago:1.00);",
file = "ex.tre", sep = "\n")
tree.primates <- read.tree("ex.tre")
X <- c(4.09434, 3.61092, 2.37024, 2.02815, -1.46968)
Y <- c(4.74493, 3.33220, 3.36730, 2.89037, 2.30259)
unlink("ex.tre") # delete the file "ex.tre"
m1 <- gls(Y ~ X, correlation=corBrownian(1, tree.primates))
corMatrix(m1$modelStruct$corStruct)
summary(m1)
m2 <- gls(Y ~ X, correlation=corMartins(1, tree.primates))
summary(m2)
corMatrix(m2$modelStruct$corStruct)
m3 <- gls(Y ~ X, correlation=corGrafen(1, tree.primates))
summary(m3)
corMatrix(m3$modelStruct$corStruct)
In [ ]:
In [ ]:
tree = Phylo.read(io.StringIO("((((Homo:0.21,Pongo:0.21):0.28,Macaca:0.49):0.13,Ateles:0.62):0.38,Galago:1.00);"),
"newick")
In [ ]:
Phylo.draw(tree)
In [ ]:
species = ["Homo","Pongo","Macaca","Ateles","Galago"]
cov = np.zeros(shape=(len(species), len(species)))
comp_matrix = np.zeros(shape=(len(species), len(species))).astype(str)
for i, spec in enumerate(species):
for j, comp in enumerate(species):
comp_matrix[i, j] = spec + ':' + comp
if i == j:
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
continue
cov[i, j] = tree.depths()[tree.common_ancestor(spec, comp)]
pd.DataFrame(cov,index=species, columns=species)
In [ ]:
X = np.array((4.09434, 3.61092, 2.37024, 2.02815, -1.46968))
Y = np.array((4.74493, 3.33220, 3.36730, 2.89037, 2.30259))
In [ ]:
p = bokeh.plotting.figure()
p.circle(X, Y)
bokeh.io.show(p)
In [ ]:
data = {
"N": len(X),
"x": X,
"y": Y,
"cov_phylo":cov,
}
sm = cmdstanpy.CmdStanModel(stan_file='spiracle_pgls_model.stan')
samples = sm.sample(data=data, sampling_iters=1000, chains=4)
samples = az.from_cmdstanpy(posterior=samples, posterior_predictive=["y_ppc"])
In [ ]:
bebi103.stan.check_all_diagnostics(samples, var_names=['b', 'a'])
In [ ]:
bokeh.io.show(
bebi103.viz.corner(
samples, pars=["a", "b"], xtick_label_orientation=np.pi / 4
)
)
In [ ]:
az.summary(samples, ['a', 'b'])
In [ ]:
np.std(np.array(samples.posterior['a']))
In [ ]:
0.7754516*np.sqrt(5)
In [ ]:
y_ppc = samples.posterior_predictive['y_ppc'].stack(
{"sample": ("chain", "draw")}
).transpose("sample", "y_ppc_dim_0")
bokeh.io.show(
bebi103.viz.predictive_regression(
y_ppc,
samples_x=X,
percentiles=[30, 50, 70, 99],
data=np.vstack((X, Y)).transpose(),
x_axis_label='ant density, k [ant/cm^2]',
y_axis_label='ant flow, q [ant/cm/s]',
x_range=[X.min()-0.1,X.max()+0.1],
)
)
In [ ]:
bebi103.stan.clean_cmdstan()
In [ ]:
cov=[[3.0, 2.0, 1.0, 1.0, 0.0],
[2.0, 3.0, 1.0, 1.0, 0.0],
[1.0, 1.0, 3.0, 2.0, 0.0],
[1.0, 1.0, 2.0, 3.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 3.0]]
In [ ]:
cov=np.array([[3.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 3.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 3.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 3.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 3.0]])*1/10
In [ ]:
x=np.array([1.02, 1.06, 0.96, 0.92, 0.89])
y=np.array([1.38, 1.41, 1.36, 1.22, 1.13])
In [ ]:
p = bokeh.plotting.figure()
p.circle(x, y)
bokeh.io.show(p)
In [ ]:
data = {
"N": len(x),
"x": x,
"y": y,
"cov_phylo":cov,
}
sm = cmdstanpy.CmdStanModel(stan_file='spiracle_pgls_model.stan')
samples = sm.sample(data=data, sampling_iters=1000, chains=4, adapt_delta=0.99)
samples = az.from_cmdstanpy(posterior=samples, posterior_predictive=["y_ppc"])
In [ ]:
bebi103.stan.check_all_diagnostics(samples, var_names=['b', 'a', 'lambda'])
In [ ]:
bokeh.io.show(
bebi103.viz.corner(
samples, pars=["a", "b", "lambda"], xtick_label_orientation=np.pi / 4
)
)
In [ ]:
az.summary(samples, ['a', 'b', 'lambda'])
In [ ]:
y_ppc = samples.posterior_predictive['y_ppc'].stack(
{"sample": ("chain", "draw")}
).transpose("sample", "y_ppc_dim_0")
bokeh.io.show(
bebi103.viz.predictive_regression(
y_ppc,
samples_x=x,
percentiles=[30, 50, 70, 99],
data=np.vstack((x, y)).transpose(),
x_axis_label='ant density, k [ant/cm^2]',
y_axis_label='ant flow, q [ant/cm/s]',
x_range=[x.min()-0.1,x.max()+0.1],
)
)
In [ ]: